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. Abstract. We briefly describe the current version of tlie PHOENIX code. 

I We tlien present some illustrative results from the modeling of Type la 

^ . and Type II supernovae, hot stars, and irradiated giant planets. Good 

fits to observations can be obtained, when account is taken for spherically 
symmetric, line-blanketed, static or expanding atmospheres. 
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1. Introduction 



PHOENIX (see Hauschildt & Baron, 1999, and references therein) is a generalized, 
stellar model atmosphere code for treating both static and moving atmospheres. 
■ The goal of PHOENIX is to be both as general as possible so that essentially all 

CD ' astrophysical objects can be modeled with a single code, and to make as few ap- 

proximations as possible. Approximations are inevitable (particularly in atomic 
data where laboratory values for most quantities are unknown); however, the 
agreement of synthetic spectra with observations across a broad class of astro- 
^! physical objects is a very good consistency check. We have modeled Planets/BDs 

P ■ (Barman et al., 2002; Allard et al., 2001; Schweitzer et al., 2001, 2002), Cool 

Stars (Hauschildt et al., 1999; Hauschildt et al., 1999), Hot Stars (/3CMa, eCMa, 
^ ; Deneb Aufdenberg et al., 1998, 1999, 2002), a-Lyra, Novae (Hauschildt et al., 

1997b; Schwarz et al., 1997), and all types of superovae (SNe la, Ib/c, IIP, lib 



Baron et al., 1995; Lentz et al., 2001a; Baron et al., 1999; Mitcheh et al., 2002). 



C3 ■ 2. The PHOENIX Code 



PHOENIX solves the radiative transfer problem by using the short-characteristic 
method (Olson et al., 1987; Olson & Kunasz, 1987) to obtain the formal solu- 
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tion of the special relativistic, spherically symmetric radiative transfer equation 
(SSRTE) along its characteristic rays. The scattering problem is solved via the 
use of a band-diagonal approximation to the discretized A-operator (Hauschildt, 
1992; Olson &: Kunasz, 1987; Hauschildt et al., 1994) as our choice of the approx- 
imate A-operator. This method can be implemented very efficiently to obtain an 
accurate solution of the SSRTE for continuum and line transfer problems using 
only modest amounts of computer resources. 

We emphasize that PHOENIX solves the radiative transfer problem including 
a full treatment of special relativistic radiative transfer in spherical geometry 
for all lines and continua. In addition we enforce the generalized condition of 
radiative equilibrium in the Lagrangian frame, including all velocity terms and 
deposition of energy from radiative decay or from external irradiation. 

We also include a full non-LTE treatment of most ions, using model atoms 
constructed from the data of Kurucz (1993, 1994a,b) and/or from the CHIANTI 
database http: //wwwsolar.nrl .navy .mil/ chianti .html. The code uses Fortran- 
95 data structures to access the different databases and model atoms from either 
or both databases can be selected at exectution time. 

Absorption and emission is treated assuming complete redistribution and 
detailed depth-dependent profiles for the lines. Fluorescence effects are included 
in the NLTE treatment. The equation of state used includes up to 26 ionization 
stages of 40 elements as well as up to 206 molecules. 

The atomic data is constructed from all relevant b-f and f-f transitions 
(Mathisen, 1984; Verner & Yakovlev, 1995; Seaton et al., 1994) as well as colli- 
sional rates obtained from laboratory measurements (where available), the Opac- 
ity Project, the Iron Project, Van Rcgemorter's formula (Van Regemorter, 1962), 
and the semi-empirical formula from Allen (1973). 

In addition to NLTE lines, lines treated in LTE and arc selected dynamically 
from the 42 x 10^ list of Kurucz (1993). A typical caknilation includes about 
2 X 10*^ lines that are treated as background opacity (Hauschildt & Baron, 1999). 
Molecular opacities (about 7 x 10^ lines) are discussed in AUard et al. (2001) 
and molecular NLTE is discussed in Schweitzer et al. (2000). 

3. The equation of radiative transfer 

The equation of transfer in spherical coordinates in the co-moving frame can be 
written (Mihalas & Mihalas, 1984): 
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h{r, jj) = 



+ 7^(1 + ^^ + 2/3/.)^ 



r 



where I{r,iJ,,X) is the specific intensity scaled by r^, r is the radial coordinate, 
// = cos 9 is the cosine of the direction angle, v is the velocity, P = v/c, 7^ = 
1/(1 — x('^) -'^) is the total extinction coefficient, % = k + Ce + Ki<f{X), and 
r/(r. A) is the emissivity. 

3.1. Numerical solution 

The details of the numerical solution are discussed in detail in Hauschildt &: 
Baron (1999) so we only sketch the basic idea here. We discretize the d/dX term 
and treat the boundary value problem for each wavelength individually, using 
an Operator splitting (OS) method. The steps are: solve along characteristics 
of the RTE, using the piecewise parabolic ansatz to calculate / for given J, and 
iterate to self-consistent solution for J using a band diagonal accelerated lambda 
operator. The eigenvalues of amplification matrix are close to unity and the use 
of operator splitting reduces the eigenvalues of the amplification matrix. 

4. Statistical Equilibrium Equations 

The radiative transfer equation depends on the opacities and emissivity, which 
in turn depend on the level populations (which in turn depend on the radia- 
tion field). Thus, we must iterate the statistical equilibrium equations, simul- 
taneously with the radiative transfer equation (the solution of the generalized 
radiative equilibrium condition is discussed in Hauschildt et al., 2003). 
The rate equations are given by (see Mihalas &; Mihalas, 1984) 



4.1. Solution of the Rate Equations 

Note that line and continuum scattering prevents the use of the A-iteration 
for the solution of the rate equations. Therefore, we use the Operator Split- 
ting method (pre-conditioning) to define a "rate operator" in analogy to the 
A-operator (Hauschildt, 1993): 



We also define an "approximate rate operator" [Rij] and write the iteration 
scheme in the form: 



j<i 




R^J = [R^J]H 



Rij = [R*j][nnew] + ([Rij] - [R*j]) [Hold]- 
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Inserting the above expression for Rij into the statistical equihbrium equa- 
tions and hnearizing we obtain, 

'^j,old[-Rji] [nnew] 



3<t 



J>1 



new J 



- rii^new I ([^^jlKld] + Cij) 



+ ^i.new ( ^ ) {[^Rji] Kid] + Cji) = 0. 



5. The Computational Problem 

The size of the computational problem is determined by: 
• input data size 

- number of atomic/ionic spectral lines: ~ 42 x 10^ 0.6 GB 

- number of diatomic molecular lines: ~ 35 x 10^ 0.5 GB 

- number of hot water vapor lines: 

- number of hot water vapor lines: 

- before 1994: w 0.035 x 10^ 

- 1994: + 6 X 10^ ^ 0.1 GB 

- 1997: + 330 X 10^ ^ 5 GB 

- 2001: + ?y 100 X 10^ ^ 1.5 GB 

- total molecular lines (May 2001): ^ 700 x 10^ 10 GB 

all lines need to be accessed in a line selection procedure which dynamically 
creates sub-lists that can be as large as the original list, but arc generally 
much smaller. This poses a significant data handling problem, which we 
have surmounted via the use of both memory and disk caches that allows 
us to trade (at run time) available memory for I/O bandwidth if required 
on a particular machine or architecture. 
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• memory /lO requirements 

— number of individual energy levels: pa 10, 000 — >^ !=a 10 MB 

— number of individual NLTE transitions: 
« 100,000 ^ PS 150 MB 

— EOS data storage « 40 MB 

— auxiliary storage ~ 50 MB 

— total memory requirement > 200 MB 

— number of individual energy levels and transitions has increased dra- 
matically — > memory requirements > 0.5 GB with the inclusion of 
more molecular species and the Chianti database. 

• (serial) CPU time 

— small for each individual point on the wavelength grid: ~ 10 . . . 100 msec 

— number of wavelength points for radiative transfer: 30,000-500,000 
(can be > 10^) 

> pa 50, 000 sec to "sweep" once through all wavelength points 

— typically, ^ 10 iterations (sweeps) are required to obtain an equilib- 
rium model 

>^6 CPU days 

— there are, literally, lOOO's of models in a typical grid . . . 

The solution is parallel computing, which we have implemented in a 
MIMD model (Hauschildt ct al., 1997a; Baron & Hauschildt, 1998; 
Hauschildt et al., 2001; Baron et al., 2003). This dramatically reduces 
wallclock time per model and makes very large scale models feasible. 

These numbers above are for models with 50 layers, with modern supercom- 
puters many models use 100 or more layers and the scaling goes approximately 
as the cube of the number of layers. 

6. Applications 

We show a few illustrative applications of PHOENIX. Figure 1 displays the ob- 
served and synthetic spectrum of the Sun and the agreement is quite good. Fig- 
ure 2 shows the observed and synthetic spectra of Vega, with excellent agreement 
in the UV and only a few features that are too strong in the optical compared 
to observations. 

Turning from stars to supernovae, we note that supernovae are fundamen- 
tally different from stars, since in stars spectra only probe to the photosphere, 
whereas a time series of supernovae spectra allow us to "peel the onion" since 
geometrical dilution due to expansion of the ejecta causes the photosphere to 
receed in mass with time. Figure 3 shows the synthetic spcctnmi of the W7 defla- 
gration model (Nomoto et al., 1984) compared to the observations of SN 1994D 
at 6 days prior to maximum light (Lentz et al., 2001a). Figure 4 shows that 
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Figure 1. The synthetic model spectrum of the Sun (thin hne) com- 
pared to observations. 



delayed detonation models well reproduce the spectra of SN 1984A — a "fast" 
SNe la (Lentz et al., 2001b; Hatano ct al., 2000). 

We have also applied PHOENIX to core collapse supernovae. Figure 5 shows 
a very good fit to the first spectrum obtained of SN 1987A in both the UV 
and optical and Figure 6 shows the PHOENIX synthetic spectrum compared to 
SN 1987A at 4.5 days after the explosion. In order to reproduce Ha, extra nickel 
mixing in the envelope was required. It is evident that the simple nickel mixing 
parameterization has destroyed the fit in the UV and clumping is probably re- 
quired to fit the entire spectrum. Figure 7 shows that excellent fits to normal 
Type IIP spectra can be obtained at early times. These fits allow us to deter- 
mine the reddening, amount of H/He mixing, amount of nickel mixing, and the 
progenitor metallicity. 



7. Conclusions 



We have described the basic design of the generalized model atmosphere code 
PHOENIX. We have shown that good, but not perfect agreement with a variety 
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Figure 2. The synthetic model spcctiiim of Vega compared to obser- 
vations (Hayes, 1985; Bohlin et al., 1990). 



of astrophysical objects can be obtained. Future work will include full time 
dependent calculations and full 3-D detailed NLTE radiative transfer. 
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Figure 3. SN 1994D on 15 March 1994 (solid line, Filippenko, 1997) 
and W7 best fit synthetic spectrum for day 14 after explosion (dashed 
line), from Lentz et al. (2001a). 
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Figure 4. Synthetic spectra for the delayed detonation model 

21c (dot-dashed lines) plotted against observed spectra for 
SN 1984A (solid lines) from Lentz et al. (2001b). 
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Figure 5. PHOENIX model spectrum (solid line) for day 1.36 (Mitchell 
et al., 2002). Important optical lines include Ha through US, and 
He I A5876. The optical spectra are taken from the CTIO archive 
(Phillips et al., 1988) and all UV spectra are from lUE (Pun et al., 
1995). 




Figure 6. PHOENIX model spectrum for SN 1987A, day 4.52, with 
gamma-ray deposition calculated assuming local deposition due to a 
constant nickel mass fraction of 1.0 x 10~^ everywhere in the envelope 
(Mitchell et al, 2002). 
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Figure 7. The observed spectrum of the Type IIP SN 1993W com- 
pared to the synthetic spectrum. The agreement is excellent. 
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